Section 4 Scenario Analysis: WCS Priority Landscapes
# turn off the s2 processing
## https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
sf::sf_use_s2(FALSE)See the USFS Wildfire Crisis Strategy (WCS) document for background information.
4.1 Scenarios Considered
There are three scenarios considered in this analysis of the constraints to mechanical forest health and fuel reduction treatments:
# mechanical mgmt allowed if:
n_sc_temp <- 3
data.frame(
scenario = 1:3
, cover = rep("Forest or Shrubland",n_sc_temp)
, protected = rep("Not Protected",n_sc_temp)
, slope = c("<40%","<60%","<60%")
, road = c("<1,000ft","<2,000ft","<2,000ft")
, riparian = c(">100ft",">100ft",">50ft")
, admin = c("No Designation","No Designation","Any Designation")
) %>%
tidyr::pivot_longer(
cols = -c(scenario)
) %>%
tidyr::pivot_wider(
names_from = scenario
, values_from = value
, names_prefix = "scenario_"
) %>%
dplyr::mutate(
name = factor(
name
, levels = c(
"cover"
, "protected"
, "slope"
, "road"
, "riparian"
, "admin"
)
, labels = c(
"NLCD Cover Type"
, "Protected or<br>IRA Status"
, "Slope"
, "Distance to<br>Nearest Road"
, "Riparian<br>Buffer"
, "Administrative<br>Designation"
)
, ordered = T
)
) %>%
dplyr::arrange(name) %>%
dplyr::mutate(
lvl = paste0("L",dplyr::row_number()-1,":")
) %>%
dplyr::relocate(lvl) %>%
# make table
kableExtra::kable(
caption = "Scenarios Considered for Determining Mechanical Treatment Operability<br>Mechanical treatment allowed if:"
, col.names = c(
""
, ""
, "Scenario 1<br>Most Constrained"
, "Scenario 2<br>Moderate"
, "Scenario 3<br>Least Constrained"
)
, escape = F
) %>%
kable_classic(full_width=T) %>%
kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE)| Scenario 1<br>Most Constrained | Scenario 2<br>Moderate | Scenario 3<br>Least Constrained | ||
|---|---|---|---|---|
| L0: | NLCD Cover Type | Forest or Shrubland | Forest or Shrubland | Forest or Shrubland |
| L1: | Protected or<br>IRA Status | Not Protected | Not Protected | Not Protected |
| L2: | Slope | <40% | <60% | <60% |
| L3: | Distance to<br>Nearest Road | <1,000ft | <2,000ft | <2,000ft |
| L4: | Riparian<br>Buffer | >100ft | >100ft | >50ft |
| L5: | Administrative<br>Designation | No Designation | No Designation | Any Designation |
The figure below illustrates the workflow used to quantify the amount of land available for mechanical forest health and risk reduction fuel treatments by considering layered operational constraints. This example illustrates scenario 2 for a 10,252 hectare (25,000 acre) fireshed project area near Rustic, Colorado in the Colorado Front Range priority landscape.

4.2 Data Preparation
4.2.1 Load WCS landscapes
WCS priority landscape spatial data from the USFS Geospatial Data Discovery site (accessed 2023-09-13).
# load data
wf_landscapes <- sf::read_sf("../data/Wildfire_Crisis_Strategy_Landscapes/Wildfire_Crisis_Strategy_Landscapes_(Feature_Layer).shp") %>%
dplyr::rename_with(tolower) %>%
sf::st_make_valid() %>%
dplyr::mutate(
hectares = (as.numeric(sf::st_area(geometry))/10000)
, Mil.Hectares = hectares %>%
scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
, acres = (as.numeric(sf::st_area(geometry))/4046.85642)
, Mil.Acres = acres %>%
scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
) %>%
dplyr::left_join(
data.frame(state = datasets::state.name, state_abb = datasets::state.abb)
, by = join_by("state")
) %>%
dplyr::mutate(
area_name = paste0(
ifelse(is.na(state_abb) | state_abb == "", "UNK", state_abb)
, ": "
, name
)
) %>%
# manual label placement
dplyr::left_join(
readr::read_csv("../data/wildfirepriority_WCS202308/area_label_placement.csv")
, by = dplyr::join_by("area_name")
)
#rename sf geom column
names(wf_landscapes)[names(wf_landscapes)==tolower(attr(wf_landscapes, "sf_column"))] = "geometry"
sf::st_geometry(wf_landscapes) = "geometry"
# set crs
transform_crs <- sf::st_crs(wf_landscapes)
# view
wf_landscapes %>% dplyr::glimpse()4.2.2 Load landscape-level constraint data
Landscape-level constraint analysis data (i.e., row unique by WCS landscape) was created via this Google Earth Engine script.
# load landscape constraint scenario data
constrained_by_scnro_ls <-
list.files("../data/wildfirepriority_WCS202308/landscapes/",pattern = "\\.csv$") %>%
purrr::keep(stringr::str_starts(.,"wfpriority_all_sc")) %>%
purrr::map(function(x){
readr::read_csv(
paste0("../data/wildfirepriority_WCS202308/landscapes/",x)
, name_repair = "universal"
, col_types = cols(.default = "c")
) %>%
dplyr::rename_with(tolower) %>%
dplyr::rename_with(make.names) %>%
dplyr::select(state,name,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>%
dplyr::mutate(scenario_id = stringr::word(x,3,sep="_") %>% readr::parse_number()) %>%
dplyr::relocate(scenario_id)
}) %>%
dplyr::bind_rows() %>%
dplyr::inner_join(
wf_landscapes %>%
sf::st_drop_geometry() %>%
dplyr::select(state,name,area_name)
, by = dplyr::join_by(state,name)
) %>%
dplyr::relocate(area_name,.after = "scenario_id") %>%
dplyr::group_by(area_name,scenario_id) %>%
dplyr::filter(dplyr::row_number()==1) %>%
dplyr::ungroup() %>%
dplyr::mutate(
dplyr::across(
tidyselect::ends_with("_m2")
, ~ as.numeric(.x) / 10000
)
, dplyr::across(
tidyselect::starts_with("pct_")
, ~ as.numeric(.x)
)
) %>%
dplyr::rename_with(
~ gsub("_m2", "_ha", .x)
, tidyselect::ends_with("_m2")
) %>%
# calculate pct reduction
dplyr::mutate(
pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
, pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
, pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
, pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
, pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
, pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
, pct_covertype_area = covertype_area_ha/feature_area_ha
# scenario description
, scenario_lab = factor(
scenario_id
, levels = 1:3
, labels = paste0("Scenario ", 1:3)
, ordered = T
) %>% forcats::fct_rev()
, scenario_desc = factor(
scenario_id
, levels = 1:3
, labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
, ordered = T
)
)4.2.3 Load fireshed project area data
The fireshed registry spatial data was obtained from the USFS Geospatial Data Discovery tool for firesheds and project areas (accessed 2023-05-03).
See this section for exploration of fireshed and project area spatial data
### fireshed spatial data
fireshed <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Fireshed/Fireshed_Registry%3A_Fireshed_(Feature_Layer).shp") %>%
sf::st_transform(transform_crs) %>%
setNames(c(
"shape_id"
, "area_ha"
, "fireshed_id"
, "fireshed_name"
, "fireshed_code"
, "fireshed_state"
, "nopas"
, "objectid"
, "fshed_id"
, "exp_total"
, "exp_usfs"
, "exp_nonfs"
, "exp_usfs_protected"
, "exp_nonfs_protected"
, "exp_usfs_managed"
, "exp_nonfs_managed"
, "exp_usfs_forest"
, "exp_nonfs_forest"
, "exp_usfs_nonforest"
, "exp_nonfs_nonforest"
, "exp_usfs_conifer"
, "exp_nonfs_conifer"
, "exp_usfs_managedforest"
, "exp_nonfs_managedforest"
, "exp_usfs_managedconifer"
, "exp_nonfs_managedconifer"
, "exp_nonfs_nonconifer_hihaz"
, "dist_vs"
, "crisis_strategy"
, "key_preformance_indicator"
, "national_usfs_rank"
, "national_all_land_rank"
, "regional_usfs_rank"
, "regional_all_land_rank"
, "start_date"
, "end_date"
, "geometry"
)) %>%
dplyr::mutate(
exposure_pct_rank = dplyr::percent_rank(exp_total)
, exposure_pct_rank_grp = dplyr::case_when(
exposure_pct_rank >= 1-0.01 ~ "Top 1%"
, exposure_pct_rank >= 1-0.05 ~ "Top 5%"
, exposure_pct_rank >= 1-0.10 ~ "Top 10%"
, exposure_pct_rank >= 1-0.25 ~ "Top 25%"
, TRUE ~ "Bottom 75%"
) %>%
factor(
levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
, ordered = T
)
# there is also a national_all_land_rank column
, ntllandrank_pct_rank = dplyr::percent_rank(-national_all_land_rank)
, ntllandrank_pct_rank_grp = dplyr::case_when(
ntllandrank_pct_rank >= 1-0.01 ~ "Top 1%"
, ntllandrank_pct_rank >= 1-0.05 ~ "Top 5%"
, ntllandrank_pct_rank >= 1-0.10 ~ "Top 10%"
, ntllandrank_pct_rank >= 1-0.25 ~ "Top 25%"
, TRUE ~ "Bottom 75%"
) %>%
factor(
levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
, ordered = T
)
, crisis_strategy = ifelse(is.na(crisis_strategy),"Not High Risk",crisis_strategy) %>%
as.factor() %>%
forcats::fct_shift()
)
#rename sf geom column
names(fireshed)[names(fireshed)==tolower(attr(fireshed, "sf_column"))] = "geometry"
sf::st_geometry(fireshed) = "geometry"
# calculate area
fireshed <- fireshed %>%
dplyr::mutate(
fireshed_area_ha = as.numeric(sf::st_area(geometry))/10000
, fireshed_area_acres = (fireshed_area_ha*10000)/4046.85642
)
## fireshed_proj_area spatial data
fireshed_proj_area <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Project_Area/Fireshed_Registry%3A_Project_Area_(Feature_Layer).shp") %>%
sf::st_transform(transform_crs) %>%
setNames(c(
"shape_id"
, "fireshed_id"
, "pa_id"
, "pa_area_ha"
, "objectid"
, "pa_id2"
, "fshed_id"
, "exp_total"
, "exp_usfs"
, "exp_nonfs"
, "exp_usfs_protected"
, "exp_nonfs_protected"
, "exp_usfs_managed"
, "exp_nonfs_managed"
, "exp_usfs_forest"
, "exp_nonfs_forest"
, "exp_usfs_nonforest"
, "exp_nonfs_nonforest"
, "exp_usfs_conifer"
, "exp_nonfs_conifer"
, "exp_usfs_managedforest"
, "exp_nonfs_managedforest"
, "exp_usfs_managedconifer"
, "exp_nonfs_managedconifer"
, "exp_nonfs_nonconifer_hihaz"
, "dist_vs"
, "pctrecentlydisturbed"
, "start_date"
, "end_date"
, "geometry"
)) %>%
dplyr::mutate(
exposure_pct_rank = dplyr::percent_rank(exp_total)
, exposure_pct_rank_grp = dplyr::case_when(
exposure_pct_rank >= 1-0.01 ~ "Top 1%"
, exposure_pct_rank >= 1-0.05 ~ "Top 5%"
, exposure_pct_rank >= 1-0.10 ~ "Top 10%"
, exposure_pct_rank >= 1-0.25 ~ "Top 25%"
, TRUE ~ "Bottom 75%"
) %>%
factor(
levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
, ordered = T
)
)
#rename sf geom column
names(fireshed_proj_area)[names(fireshed_proj_area)==tolower(attr(fireshed_proj_area, "sf_column"))] = "geometry"
sf::st_geometry(fireshed_proj_area) = "geometry"
# calculate area
fireshed_proj_area <- fireshed_proj_area %>%
dplyr::mutate(
pa_area_ha = as.numeric(sf::st_area(geometry))/10000
, pa_area_acres = (pa_area_ha*10000)/4046.85642
) %>%
# JOIN WITH FIRESHED DATA
dplyr::inner_join(
fireshed %>%
sf::st_drop_geometry() %>%
dplyr::select(fireshed_id, crisis_strategy, exp_total
, exposure_pct_rank, exposure_pct_rank_grp
) %>%
dplyr::rename(exposure_total=exp_total) %>%
dplyr::rename_with(
~ paste0("fireshed_",.x)
, -c(fireshed_id)
)
, by = dplyr::join_by(fireshed_id)
) %>%
dplyr::select(pa_id,pa_area_ha
,exp_total,exposure_pct_rank,exposure_pct_rank_grp
, tidyselect::starts_with("fireshed_")
) %>%
dplyr::rename(exposure_total=exp_total) %>%
dplyr::rename_with(
~ paste0("pa_", .x, recycle0 = TRUE)
, tidyselect::starts_with("exp")
)Fireshed project area constraint analysis data (i.e., row unique by WCS landscape + FPA) was created via this Google Earth Engine script
# load fireshed constraint scenario data
constrained_by_scnro_ls_pa <-
list.files("../data/wildfirepriority_WCS202308/fireshed/",pattern = "\\.csv$") %>%
purrr::keep(stringr::str_starts(.,"wfpriority_fireshed_sc")) %>%
purrr::map(function(x){
readr::read_csv(
paste0("../data/wildfirepriority_WCS202308/fireshed/",x)
, name_repair = "universal"
, col_types = cols(.default = "c")
) %>%
dplyr::rename_with(tolower) %>%
dplyr::rename_with(make.names) %>%
dplyr::select(state,name,pa_id,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>%
dplyr::mutate(scenario_id = stringr::word(x,3,sep="_") %>% readr::parse_number()) %>%
dplyr::relocate(scenario_id)
}) %>%
dplyr::bind_rows() %>%
dplyr::inner_join(
wf_landscapes %>%
sf::st_drop_geometry() %>%
dplyr::select(state,name,area_name)
, by = dplyr::join_by(state,name)
) %>%
dplyr::relocate(area_name,.after = "scenario_id") %>%
dplyr::mutate(
dplyr::across(
tidyselect::ends_with("_m2")
, ~ as.numeric(.x) / 10000
)
, dplyr::across(
tidyselect::starts_with("pct_")
, ~ as.numeric(.x)
)
) %>%
dplyr::rename_with(
~ gsub("_m2", "_ha", .x)
, tidyselect::ends_with("_m2")
) %>%
dplyr::group_by(area_name,pa_id,scenario_id) %>%
dplyr::arrange(desc(pct_pa_intrsct)) %>%
dplyr::filter(dplyr::row_number()==1) %>%
dplyr::ungroup() %>%
dplyr::filter(pct_pa_intrsct>=0.00001) %>%
# calculate pct reduction
dplyr::mutate(
pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
, pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
, pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
, pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
, pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
, pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
, pct_covertype_area = covertype_area_ha/feature_area_ha
# calculate level of constraint
, cnstrnt_lvl = dplyr::case_when(
-pct_rdctn_total > 0.8 ~ 1
, -pct_rdctn_total >= 0.6 ~ 2
, -pct_rdctn_total >= 0.0 ~ 3
)
, cnstrnt_class = factor(
cnstrnt_lvl
, levels = 1:3
, labels = c("high constraint", "med. constraint", "low constraint")
, ordered = T
) %>% forcats::fct_rev()
, rmn_cnstrnt_class = factor(
cnstrnt_lvl
, levels = 1:3
, labels = c("0–19% treatable", "20–40% treatable", ">40% treatable")
, ordered = T
) %>% forcats::fct_rev()
# scenario description
, scenario_lab = factor(
scenario_id
, levels = 1:3
, labels = paste0("Scenario ", 1:3)
, ordered = T
) %>% forcats::fct_rev()
, scenario_desc = factor(
scenario_id
, levels = 1:3
, labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
, ordered = T
)
) %>%
# join with fireshed data
dplyr::inner_join(
fireshed_proj_area %>%
sf::st_drop_geometry() %>%
dplyr::select(pa_id, tidyselect::starts_with("pa_exposure"), tidyselect::starts_with("fireshed_")) %>%
dplyr::mutate(pa_id=as.character(pa_id))
, by = dplyr::join_by(pa_id)
)4.3 Geographic Delineation Description
Timeline of the research utilized to motivate and support the creation of the WCS:
- 2018 August - Toward share stewardship across the landscapes: an outcome-based investment strategy created a hotspot map of the source of wildfire risk to communities showing that 80% of community exposure from NFS lands originates on 20% of the NFS area.
- 2019 February - Evers et al. analyzed community wildfire exposure from national forests by associating conditions that affect exposure in the areas where wildfires ignite to conditions where exposure likely occurs
- 2019 May - Ager et al. mapped community firesheds by creating a continuous smoothed surface of predicted structure exposure from all FSim ignitions (see also Palaiologou et al. 2019 and Ager et al. 2019)
- 2020 - Evers et al. release Fireshed Registry as a geospatial dashboard and decision tool
- 2021 May - Ager et al. 2021 official release of the Fireshed Registry as a geospatial information system that organizes landscape risk to developed areas into a planning framework (see figure below)
- 2022 January - US Department of Agriculture, Forest Service launches a 10-year Wildfire Crisis Strategy
- 2022 April - 10 initial landscapes selected for WCS prioritization
- 2023 January - additional 11 landscapes selected for WCS prioritization

Nested spatial framework for firesheds. Each scale has specific functionality in terms of the planning processes. Firesheds are the broad scale unit of prioritization, but project areas within them are also prioritized as part of implementation of treatments. Project areas are roughly the size that national forests use for conducting vegetation and fuel management projects. The relative variation among firesheds compared to variation within them controls the relative emphasis on selecting firesheds versus individual planning areas. Ager et al. 2021
#define basemap
states_temp <- USAboundaries::us_states(
states = c(
# usfs region 1-6 states
"MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
# , "KS","NE","SD","ND"
)
) %>%
sf::st_transform(transform_crs)
cities_temp <- USAboundaries::us_cities(
states = c(
# usfs region 1-6 states
"MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
# , "KS","NE","SD","ND"
)
) %>%
sf::st_transform(transform_crs) %>%
dplyr::filter(
toupper(city) %in% c(
"DENVER"
, "TUSCON"
, "SALT LAKE CITY"
, "LAS VEGAS"
, "SAN DIEGO"
, "LOS ANGELES"
, "FRESNO"
, "SAN FRANCISCO"
, "SACRAMENTO"
, "PORTLAND"
, "SEATTLE"
, "ALBUQUERQUE"
, "RENO"
, "BOISE"
, "HELENA"
, "BILLINGS"
)
| (toupper(city) == "PHOENIX" & state_abbr == "AZ")
)
# plot basemap
plt_basemap <-
ggplot() +
geom_sf(
data = states_temp
, fill = NA
, color = "gray66"
) +
geom_sf_text(
data = states_temp
, mapping = aes(label = stusps)
, color = "gray66"
, size = 3
) +
geom_sf(
data = cities_temp
, shape = 1
, size = 1
, color = "gray44"
) +
geom_sf_text(
data = cities_temp
, mapping = aes(label = city)
, color = "gray44"
, size = 2
# , hjust = -0.15
, hjust = 0.5
, vjust = -0.25
) +
coord_sf(expand = F) +
theme_void()
# plot basemap simple
plt_basemap_simple <-
ggplot() +
geom_sf(
data = states_temp
, fill = NA
, color = "gray66"
) +
coord_sf(expand = F) +
theme_void()4.3.1 Fireshed risk to communities (2021)
ggplot() +
geom_sf(
data = fireshed
, mapping = aes(fill=ntllandrank_pct_rank_grp)
) +
geom_sf(
data = USAboundaries::us_states() %>%
dplyr::filter(!stusps %in% c("AK","HI","PR")) %>%
sf::st_transform(transform_crs)
, fill = NA
, color = "gray88"
) +
scale_fill_manual(values=c("firebrick","orange3","khaki","seagreen3","royalblue")) +
coord_sf(expand = F) +
labs(
fill = "Fireshed exposure ranks"
) +
theme_void() +
theme(
legend.position = "top"
, legend.direction = "horizontal"
, legend.title = element_text(size = 8)
, legend.text = element_text(size = 7)
)
National map of the 7,688 firesheds created from community wildfire transmission data (Evers et al. 2020). The fireshed boundaries were created with a process that delineates hotspots of fire transmission to buildings in adjacent or nearby communities. See the Methods section for details on delineating firesheds. Figure 2 from Ager et al. 2021 (p. 7)
4.3.2 Wilfire Crisis Strategy 21 Priority Landscapes (2022-2023)
plt_fshed <-
plt_basemap +
geom_sf(
data = fireshed %>%
sf::st_filter(
wf_landscapes %>%
dplyr::select(area_name)
, .predicate=st_intersects
)
, aes(fill=crisis_strategy)
, alpha = 0.8
) +
geom_sf(
data = wf_landscapes
, mapping = aes(
color = paste0(investment, " investment")
# , linetype = paste0(investment, " investment")
)
, fill=NA
# , color="black"
, lwd=0.7
) +
geom_text(
data = wf_landscapes %>%
sf::st_drop_geometry() %>%
sf::st_as_sf(coords = c("lon","lat")) %>%
sf::st_set_crs(4326) %>%
sf::st_transform(transform_crs)
, aes(
label = dplyr::case_when(
stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
, T ~
stringr::str_wrap(
dplyr::case_when(
stringr::str_count(name, "\\w+")<2
~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
, TRUE ~ stringr::word(name,1,2)
)
, 7
)
)
, color = paste0(investment, " investment")
, geometry = geometry
, fontface = "bold.italic"
)
, stat = "sf_coordinates"
, size = 2.5
# , seed = 11
# , min.segment.length = 0
, show.legend = F
) +
scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
scale_color_manual(values = c("gray11", "navy")) + # midnightblue
labs(
fill = "Fireshed\nHigh-Risk Status"
, color = "Landscapes"
# , linetype = "Landscapes"
) +
# ggspatial::annotation_scale(location = "br", style = "ticks", pad_y = unit(0.7, "cm"), width_hint = 0.05) +
theme(
legend.position = c(0.87,0.87) # "top"
# , legend.direction = "horizontal"
, legend.title = element_text(size = 8, face = "bold")
, legend.text = element_text(size = 8)
) +
guides(
color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
, fill = guide_legend(order = 2, override.aes = list(lwd = 0))
)
# plot
plt_fshed
ggplot2::ggsave(
filename = paste0("../data/plt_map_ls_fshed.jpeg")
, plot = ggplot2::last_plot()
, width = 8.5
, height = 8.5
, units = "in"
, dpi = "print"
, bg = "white"
)The boundaries of the 21 priority landscapes roughly follow the boundaries of “firesheds” prioritized to reduce wildfire transmission to developed areas (Figure 1). Firesheds are geographic delineations with an average area of approximately 250,000 acres (~100,000 hectares) that were created to organize the landscape into units for managing wildfire risk to communities.
4.3.3 Wilfire Crisis Strategy Treatment Objective
plt_basemap +
geom_sf(
data = fireshed %>%
dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>%
sf::st_filter(
wf_landscapes %>%
dplyr::select(area_name)
, .predicate=st_intersects
)
, aes(fill=crisis_strategy)
, alpha = 0.8
) +
geom_sf(
data = wf_landscapes
, mapping = aes(color = paste0(investment, " investment"))
, fill=NA
# , color="black"
, lwd=0.7
)+
geom_text(
data = wf_landscapes %>%
sf::st_drop_geometry() %>%
sf::st_as_sf(coords = c("lon","lat")) %>%
sf::st_set_crs(4326) %>%
sf::st_transform(transform_crs)
, aes(
label = dplyr::case_when(
stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
, T ~
stringr::str_wrap(
dplyr::case_when(
stringr::str_count(name, "\\w+")<2
~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
, TRUE ~ stringr::word(name,1,2)
)
, 7
)
)
, color = paste0(investment, " investment")
, geometry = geometry
, fontface = "bold.italic"
)
, stat = "sf_coordinates"
, size = 2.5
# , seed = 11
# , min.segment.length = 0
, show.legend = F
) +
scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
scale_color_manual(values = c("black", "navy")) +
# ggspatial::annotation_scale(location = "br", style = "ticks", pad_y = unit(0.7, "cm"), width_hint = 0.05) +
labs(
fill = "Fireshed\nHigh-Risk Status"
, color = "Landscapes"
) +
theme(
legend.position = c(0.87,0.87) # "top"
# , legend.direction = "horizontal"
, legend.title = element_text(size = 8, face = "bold")
, legend.text = element_text(size = 8)
) +
guides(
color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
, fill = guide_legend(order = 2, override.aes = list(lwd = 0))
)
Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).
4.3.4 Fireshed Project Areas
name_temp <- "Enchanted Circle"
# pa map
plt_fshed_pa <-
ggplot() +
geom_sf(data =
fireshed_proj_area %>%
dplyr::mutate(pa_id=as.character(pa_id)) %>%
dplyr::inner_join(
constrained_by_scnro_ls_pa %>%
dplyr::filter(
# fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
name == name_temp
) %>%
dplyr::select(pa_id)
, by = dplyr::join_by(pa_id)
, multiple = "all"
)
, aes(fill = fireshed_crisis_strategy, color = "Project Area\nboundary")
, alpha = 0.8
, lwd = 1.2
) +
geom_sf(
data = fireshed %>%
# dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>%
sf::st_filter(
wf_landscapes %>%
dplyr::filter(name == name_temp) %>%
dplyr::select(area_name)
, .predicate=st_intersects
)
, aes(color = "Fireshed\nboundary")
, fill = NA
, lwd = 1.2
, linetype = "dashed"
) +
geom_sf(
data = wf_landscapes %>% dplyr::filter(name == name_temp)
, mapping = aes(color = "WCS Landscape\nboundary")
, fill=NA
, lwd = 0.9
)+
scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
scale_color_manual(values = c("gray44","skyblue2", "black")) +
labs(
fill = "Fireshed\nHigh-Risk Status"
, color = "Spatial\nFramework"
, subtitle = name_temp
) +
coord_sf(expand = F) +
ggspatial::annotation_scale(location = "br", style = "ticks", pad_x = unit(0.1, "cm"), pad_y = unit(0.25, "cm")) +
theme_void() +
theme(
legend.position = c(1.08,0.79) # "top"
# , legend.direction = "horizontal"
, legend.title = element_text(size = 8, face = "bold")
, legend.text = element_text(size = 10)
, plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
) +
guides(
color = guide_legend(order = 1, override.aes = list(size = 2.7, linewidth = 6, fill = NA, linetype = 1)) #, linetype = c(5,1,1)))
, fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
)# pa map w/in
plt_fshed_pa_incl <- ggplot() +
geom_sf(data =
fireshed_proj_area %>%
dplyr::mutate(pa_id=as.character(pa_id)) %>%
dplyr::inner_join(
constrained_by_scnro_ls_pa %>%
dplyr::filter(
# fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
name == name_temp
) %>%
dplyr::mutate(
# is_25in = ifelse(pct_pa_intrsct>=0.25,"≥25%","<25%")
is_25in = ifelse(pct_pa_intrsct>=0.25,"included","excluded")
) %>%
dplyr::select(
pa_id
, is_25in
)
, by = dplyr::join_by(pa_id)
, multiple = "all"
)
, aes(fill = is_25in, color = "Project Area\nboundary")
, alpha = 0.9
, lwd = 1.2
) +
geom_sf(
data = fireshed %>%
# dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>%
sf::st_filter(
wf_landscapes %>%
dplyr::filter(name == name_temp) %>%
dplyr::select(area_name)
, .predicate=st_intersects
)
, aes(color = "Fireshed\nboundary")
, fill = NA
, lwd = 1.2
, linetype = "dashed"
) +
geom_sf(
data = wf_landscapes %>% dplyr::filter(name == name_temp)
, mapping = aes(color = "WCS Landscape\nboundary")
, fill=NA
, lwd = 0.9
)+
scale_fill_manual(values = viridis::viridis(n=4, direction = -1, alpha = 0.9)[c(1,3)]) +
scale_color_manual(values = c("gray44","skyblue2", "black")) +
labs(
fill = "Project Area\nInclusions"
# fill = "Included in this analysis\nif ≥25% area w/in landscape"
, color = "Spatial\nFramework"
, subtitle = name_temp
) +
coord_sf(expand = F) +
ggspatial::annotation_scale(location = "br", style = "ticks", pad_x = unit(0.1, "cm"), pad_y = unit(0.25, "cm")) +
theme_void() +
theme(
legend.position = c(1.08,0.825) # "top"
# , legend.direction = "horizontal"
, legend.title = element_text(size = 8, face = "bold")
, legend.text = element_text(size = 10)
, plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
) +
guides(
color = guide_legend(order = 1, override.aes = list(size = 2.7, linewidth = 6, fill = NA, linetype = 1)) #, linetype = c(5,1,1)))
, fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
)
# combine
cowplot::plot_grid(
plotlist = list(
NULL
, plt_fshed_pa +
theme(
plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
, plot.margin = margin(0.5,3,0.5,0.5, "cm")
)
, NULL
, plt_fshed_pa_incl +
theme(
plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
, plot.margin = margin(0.5,3,0.5,0.5, "cm")
)
, NULL
)
, nrow = 1
, rel_widths = c(0.01,1,0.01,1,0.01)
, labels = c("","A","","B","")
, label_size = 17
, label_fontface = "bold"
, label_colour = "gray11"
# , label_y = 0.99
) +
theme(plot.background = element_rect(fill = "white", color = NA))
# export
ggplot2::ggsave(
filename = paste0("../data/plt_ex_spatial_strctr.jpeg")
, plot = ggplot2::last_plot()
, width = 11
, height = 5.8
, units = "in"
, dpi = "print"
)Nested within firesheds are project areas of approximately 10,000 ha in size which represent the geographic unit at which vegetation and fuel management projects are planned (Ager et al. 2021).
Fireshed project areas that did not have at least 25% of area within the boundary of the Strategy priority landscape were excluded from this analysis. This cutoff was implemented under the assumption that with less than 25% of area available, treatment alone would not substantially affect wildfire behavior across the fireshed (North et al. 2015). For each fireshed project area included in this analysis, the entire area within the project area boundary was used to calculate the percent of mechanically available acreage including area outside of the priority landscape boundary.
4.4 Landscape-level analysis
# filter fireshed project area
constrained_by_scnro_ls_pa <- constrained_by_scnro_ls_pa %>%
dplyr::filter(pct_pa_intrsct>=0.25)4.4.1 Reduction Treatable Area Table
# scnum_temp<-2
tbl_temp <- constrained_by_scnro_ls %>%
# dplyr::filter(scenario_id==scnum_temp) %>%
dplyr::mutate(
dplyr::across(
tidyselect::ends_with("_ha")
, ~ scales::comma(.x, suffix = "k", scale = 1e-3, accuracy = 1)
)
, dplyr::across(
tidyselect::starts_with("pct_")
, ~ scales::percent(.x, accuracy = .1)
)
) %>%
dplyr::select(
area_name
, scenario_lab
, covertype_area_ha
, pct_covertype_area
, pct_rdctn1_protected
, pct_rdctn2_slope
, pct_rdctn3_roads
, pct_rdctn4_riparian
, pct_rdctn5_administrative
, rmn5_administrative_area_ha
, pct_rmn5_administrative
) %>%
dplyr::arrange(area_name,desc(scenario_lab))
# make table
kableExtra::kable(
tbl_temp %>% dplyr::select(-c(area_name))
# , caption = paste0("<b><font color=navy>Scenario "
# ,scnum_temp
# ,"</font></b><br>percent reduction of different types of constraints on mechanical treatment"
# )
, caption = "<b>Wildfire Crisis Strategy 21 Priority Landscapes</b><br>percent reduction of different types of constraints on mechanical treatment"
, col.names = c(
""
, "Forest+Shrub\n(ha)", "Forest+Shrub\n(%)"
, "Protected"
, "Slope\nSteepness"
, "Road\nDistance"
, "Riparian\nBuffer"
, "Administrative"
, "Mechanically\nAvailable (ha)"
, "Mechanically\nAvailable (%)"
)
, escape = F
) %>%
add_header_above(c(" " = 1, "Area of\nPriority Landscape"=2, "Constraint\nLeast Flexible to Most Flexible" = 5, " " = 2)) %>%
kable_classic(full_width=T) %>%
pack_rows(index = table(forcats::fct_inorder(tbl_temp$area_name))) %>%
kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%
kableExtra::scroll_box(width = "740px")| Forest+Shrub (ha) | |Forest+Shrub (% | |Protecte | |Slope Steepne | s |Road Dista | ce |Riparian Bu | fer |Administra | ive |Mechanically Available | (ha) |Mechanically Availab | |
|---|---|---|---|---|---|---|---|---|---|
| AZ: 4FRI | |||||||||
| Scenario 1 | 2,128k | 90.4% | -12.4% | -7.5% | -26.7% | -10.2% | -15.5% | 589k | 27.7% |
| Scenario 2 | 2,128k | 90.4% | -12.4% | -2.1% | -13.3% | -14.0% | -19.9% | 813k | 38.2% |
| Scenario 3 | 2,128k | 90.4% | -12.4% | -2.1% | -13.3% | -9.9% | 0.0% | 1,325k | 62.3% |
| AZ: Prescott | |||||||||
| Scenario 1 | 275k | 87.2% | -12.9% | -16.1% | -25.1% | -11.9% | -4.2% | 82k | 29.8% |
| Scenario 2 | 275k | 87.2% | -12.9% | -3.3% | -13.2% | -18.6% | -6.8% | 124k | 45.2% |
| Scenario 3 | 275k | 87.2% | -12.9% | -3.3% | -13.2% | -13.2% | 0.0% | 158k | 57.4% |
| AZ: San Carlos Apache Tribal Forest Protection | |||||||||
| Scenario 1 | 1,119k | 91.7% | -9.3% | -19.7% | -49.8% | -5.8% | -1.1% | 161k | 14.4% |
| Scenario 2 | 1,119k | 91.7% | -9.3% | -5.0% | -46.6% | -10.1% | -2.2% | 300k | 26.8% |
| Scenario 3 | 1,119k | 91.7% | -9.3% | -5.0% | -46.6% | -7.3% | 0.0% | 357k | 31.9% |
| CA: North Yuba | |||||||||
| Scenario 1 | 139k | 96.5% | -11.6% | -33.4% | -7.8% | -14.7% | -2.3% | 42k | 30.2% |
| Scenario 2 | 139k | 96.5% | -11.6% | -13.9% | -3.2% | -23.0% | -3.3% | 63k | 45.0% |
| Scenario 3 | 139k | 96.5% | -11.6% | -13.9% | -3.2% | -16.7% | 0.0% | 76k | 54.6% |
| CA: Plumas Community Protection | |||||||||
| Scenario 1 | 101k | 93.3% | 0.0% | -23.2% | -19.8% | -21.5% | -1.9% | 34k | 33.6% |
| Scenario 2 | 101k | 93.3% | 0.0% | -5.2% | -8.1% | -33.9% | -2.9% | 50k | 49.9% |
| Scenario 3 | 101k | 93.3% | 0.0% | -5.2% | -8.1% | -24.4% | 0.0% | 63k | 62.3% |
| CA: Southern California Fireshed Risk Reduction Strategy | |||||||||
| Scenario 1 | 1,408k | 86.0% | -57.0% | -21.0% | -8.3% | -2.3% | -1.4% | 139k | 9.9% |
| Scenario 2 | 1,408k | 86.0% | -57.0% | -9.6% | -6.5% | -4.2% | -3.0% | 278k | 19.7% |
| Scenario 3 | 1,408k | 86.0% | -57.0% | -9.6% | -6.5% | -3.0% | 0.0% | 336k | 23.9% |
| CA: Stanislaus | |||||||||
| Scenario 1 | 116k | 94.3% | -2.9% | -26.7% | -10.6% | -24.7% | -0.6% | 40k | 34.5% |
| Scenario 2 | 116k | 94.3% | -2.9% | -10.7% | -3.9% | -35.8% | -0.8% | 53k | 46.0% |
| Scenario 3 | 116k | 94.3% | -2.9% | -10.7% | -3.9% | -26.1% | 0.0% | 66k | 56.5% |
| CA: Trinity Forest Health and Fire-Resilient Rural Communities | |||||||||
| Scenario 1 | 554k | 83.5% | -21.4% | -37.8% | -10.4% | -8.7% | -9.2% | 69k | 12.4% |
| Scenario 2 | 554k | 83.5% | -21.4% | -12.4% | -7.8% | -17.3% | -18.5% | 126k | 22.7% |
| Scenario 3 | 554k | 83.5% | -21.4% | -12.4% | -7.8% | -12.5% | 0.0% | 254k | 45.9% |
| CO: Colorado Front Range | |||||||||
| Scenario 1 | 1,052k | 72.6% | -20.4% | -19.9% | -20.0% | -11.3% | -1.6% | 281k | 26.7% |
| Scenario 2 | 1,052k | 72.6% | -20.4% | -5.6% | -12.0% | -17.3% | -3.2% | 436k | 41.4% |
| Scenario 3 | 1,052k | 72.6% | -20.4% | -5.6% | -12.0% | -12.3% | 0.0% | 522k | 49.6% |
| ID: Nez Perce-Clearwater-Lower Salmon | |||||||||
| Scenario 1 | 700k | 89.0% | -43.0% | -19.7% | -9.1% | -1.4% | -0.3% | 185k | 26.5% |
| Scenario 2 | 700k | 89.0% | -43.0% | -6.7% | -5.5% | -2.3% | -0.7% | 293k | 41.8% |
| Scenario 3 | 700k | 89.0% | -43.0% | -6.7% | -5.5% | -1.6% | 0.0% | 302k | 43.2% |
| ID: Southwest Idaho | |||||||||
| Scenario 1 | 611k | 87.6% | -13.3% | -28.2% | -10.4% | -2.0% | -0.3% | 279k | 45.8% |
| Scenario 2 | 611k | 87.6% | -13.3% | -7.0% | -6.6% | -2.8% | -0.8% | 425k | 69.5% |
| Scenario 3 | 611k | 87.6% | -13.3% | -7.0% | -6.6% | -2.0% | 0.0% | 435k | 71.2% |
| MT: Kootenai Complex | |||||||||
| Scenario 1 | 372k | 91.4% | -21.4% | -15.0% | -7.8% | -6.9% | -10.6% | 142k | 38.3% |
| Scenario 2 | 372k | 91.4% | -21.4% | -3.6% | -2.1% | -8.6% | -13.8% | 188k | 50.6% |
| Scenario 3 | 372k | 91.4% | -21.4% | -3.6% | -2.1% | -6.1% | 0.0% | 249k | 66.9% |
| NM: Enchanted Circle | |||||||||
| Scenario 1 | 506k | 85.6% | -8.5% | -21.5% | -33.4% | -7.5% | -2.2% | 136k | 26.8% |
| Scenario 2 | 506k | 85.6% | -8.5% | -5.2% | -25.7% | -11.7% | -3.7% | 229k | 45.2% |
| Scenario 3 | 506k | 85.6% | -8.5% | -5.2% | -25.7% | -8.2% | 0.0% | 265k | 52.3% |
| NV: Sierra and Elko Fronts | |||||||||
| Scenario 1 | 1,001k | 73.1% | -26.6% | -9.8% | -26.8% | -6.2% | -0.9% | 297k | 29.6% |
| Scenario 2 | 1,001k | 73.1% | -26.6% | -2.0% | -15.7% | -8.9% | -1.5% | 453k | 45.3% |
| Scenario 3 | 1,001k | 73.1% | -26.6% | -2.0% | -15.7% | -6.3% | 0.0% | 495k | 49.4% |
| OR: Central Oregon | |||||||||
| Scenario 1 | 959k | 87.1% | -7.0% | -2.8% | -15.3% | -5.4% | -10.0% | 571k | 59.5% |
| Scenario 2 | 959k | 87.1% | -7.0% | -0.6% | -5.0% | -6.4% | -11.3% | 669k | 69.7% |
| Scenario 3 | 959k | 87.1% | -7.0% | -0.6% | -5.0% | -4.5% | 0.0% | 796k | 83.0% |
| OR: Klamath River Basin | |||||||||
| Scenario 1 | 2,675k | 77.0% | -19.3% | -14.9% | -17.7% | -5.4% | -4.4% | 1,023k | 38.3% |
| Scenario 2 | 2,675k | 77.0% | -19.3% | -4.9% | -7.9% | -8.3% | -7.8% | 1,386k | 51.8% |
| Scenario 3 | 2,675k | 77.0% | -19.3% | -4.9% | -7.9% | -5.9% | 0.0% | 1,659k | 62.0% |
| OR: Mount Hood Forest Health and Fire-Resilient Communities | |||||||||
| Scenario 1 | 326k | 76.3% | -20.1% | -14.3% | -17.1% | -9.0% | -15.7% | 78k | 23.8% |
| Scenario 2 | 326k | 76.3% | -20.1% | -4.1% | -7.5% | -13.2% | -20.1% | 114k | 35.0% |
| Scenario 3 | 326k | 76.3% | -20.1% | -4.1% | -7.5% | -9.5% | 0.0% | 191k | 58.8% |
| UT: Pine Valley | |||||||||
| Scenario 1 | 153k | 94.2% | -37.2% | -6.9% | -21.5% | -6.7% | 0.0% | 43k | 27.8% |
| Scenario 2 | 153k | 94.2% | -37.2% | -1.4% | -10.5% | -9.2% | 0.0% | 64k | 41.7% |
| Scenario 3 | 153k | 94.2% | -37.2% | -1.4% | -10.5% | -6.5% | 0.0% | 68k | 44.4% |
| UT: Wasatch | |||||||||
| Scenario 1 | 381k | 89.5% | -43.0% | -12.1% | -14.8% | -4.9% | -0.3% | 95k | 25.0% |
| Scenario 2 | 381k | 89.5% | -43.0% | -3.7% | -7.9% | -6.7% | -0.5% | 146k | 38.3% |
| Scenario 3 | 381k | 89.5% | -43.0% | -3.7% | -7.9% | -4.8% | 0.0% | 155k | 40.7% |
| WA: Central Washington Initiative | |||||||||
| Scenario 1 | 893k | 69.8% | -25.0% | -29.5% | -11.5% | -8.1% | -12.7% | 117k | 13.1% |
| Scenario 2 | 893k | 69.8% | -25.0% | -9.5% | -7.0% | -13.4% | -22.2% | 204k | 22.9% |
| Scenario 3 | 893k | 69.8% | -25.0% | -9.5% | -7.0% | -9.7% | 0.0% | 436k | 48.8% |
| WA: Colville Northeast Washington Vision | |||||||||
| Scenario 1 | 656k | 92.8% | -20.1% | -18.5% | -16.0% | -10.1% | -1.2% | 223k | 34.1% |
| Scenario 2 | 656k | 92.8% | -20.1% | -3.8% | -7.5% | -14.4% | -2.0% | 342k | 52.2% |
| Scenario 3 | 656k | 92.8% | -20.1% | -3.8% | -7.5% | -10.4% | 0.0% | 382k | 58.3% |
4.4.2 Change in Available by Scenario
constrained_by_scnro_ls %>%
dplyr::mutate(scenario_lab = scenario_lab %>% forcats::fct_rev()) %>%
ggplot(
mapping = aes(x=scenario_lab,y=pct_rmn5_administrative,group=area_name)
) +
geom_line(mapping=aes(color = area_name), linewidth = 1.5) +
geom_label(
mapping=aes(label = scales::percent(pct_rmn5_administrative, scale = 100, accuracy = 1))
, color = "black"
, size = 3
, label.padding = unit(0.15, "lines")
) +
facet_wrap(facets = vars(area_name)
, ncol = 3
, labeller = label_wrap_gen(width = 35, multi_line = TRUE)
, scales = "free_x"
) +
scale_color_viridis_d(option = "turbo", alpha = 0.8) +
scale_y_continuous(expand = expansion(mult = c(0.1, 0.1)), labels = scales::percent_format(), breaks = scales::extended_breaks(6)) +
# scale_x_discrete(labels = scales::label_wrap(20)) +
labs(
x = "Most Restrictive \U2192 Least Restrictive"
, y = "Percent Treatable Area Remaining"
) +
theme_light() +
theme(
legend.position = "none"
, axis.text.y = element_text(size=7)
, axis.text.x = element_text(size=7.5)
, strip.text = element_text(color = "black", face = "bold", size = 10)
)
4.4.3 Reduction by Constraint
# reshape
constrained_by_scnro_ls_long <- constrained_by_scnro_ls %>%
dplyr::select(state, name, area_name, tidyselect::starts_with("scenario_"), tidyselect::starts_with("pct_rdctn")) %>%
tidyr::pivot_longer(
cols = tidyselect::starts_with("pct_rdctn")
, names_to = "constraint"
, values_to = "pct_rdctn"
, names_prefix = "pct_rdctn"
, values_drop_na = F
) %>%
tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>%
dplyr::mutate(
constraint_lvl = as.numeric(constraint_lvl)
, constraint = factor(
constraint
, ordered = TRUE
, levels = c(
"protected"
, "slope"
, "roads"
, "riparian"
, "administrative"
, "total"
)
, labels = c(
"Protected"
, "Slope\nSteepness"
, "Road\nDistance"
, "Riparian\nBuffer"
, "Administrative"
, "Total"
)
) %>% forcats::fct_rev()
) %>%
dplyr::left_join(
constrained_by_scnro_ls %>% dplyr::select(state,name,scenario_id,pct_rdctn_total)
, by = join_by(state,name,scenario_id)
)# color pallete removing dark blue
n_temp <- ((constrained_by_scnro_ls_long$constraint %>% unique() %>% length())-1)*2
plasma_custom <- viridis::plasma(n = n_temp)[seq(2,n_temp,2)]
# scales::show_col(plasma_custom)
# plasma_custom
# plot
ggplot(
data = constrained_by_scnro_ls_long %>%
dplyr::filter(constraint!="Total") %>%
dplyr::mutate(
dplyr::across(
c(scenario_lab, area_name)
, ~ forcats::fct_rev(.x)
)
)
, mapping = aes(y = area_name, x = pct_rdctn)
) +
geom_col(
mapping = aes(fill = constraint)
, color = NA, width = 0.7
, alpha = 0.7
) +
geom_text(
mapping = aes(
label = scales::percent(ifelse(pct_rdctn<0.11*-1,pct_rdctn,NA), accuracy = 1)
, group = constraint
, fontface = "bold"
)
, position = position_stack(vjust = 0.5)
, color = "black", size = 2.3
) +
geom_text(
data = constrained_by_scnro_ls_long %>%
dplyr::filter(constraint=="Total") %>%
dplyr::mutate(
dplyr::across(
c(scenario_lab, area_name)
, ~ forcats::fct_rev(.x)
)
)
, mapping = aes(
y = area_name, x = pct_rdctn_total
, label = scales::percent(pct_rdctn_total, accuracy = 1)
)
, color = "black", size = 3.5
, hjust = -0.1
) +
facet_grid(cols = vars(scenario_lab)) +
# scale_fill_viridis_d(option = "plasma") +
scale_fill_manual(values = plasma_custom) +
scale_x_reverse(expand = expansion(mult = c(0.04, .3)),labels = scales::percent_format()) +
labs(
fill = ""
, y = ""
, x = "Constraint Reduction in Treatable Area (%)"
) +
theme_light() +
theme(
legend.position = "top"
, legend.title = element_text(size=7)
, axis.title = element_text(size=9)
, axis.text.x = element_blank() #
, strip.text = element_text(color = "black", face = "bold")
) +
guides(
fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
)
ggplot2::ggsave(
filename = paste0("../data/plt_cnstrnt_rdctn_ls.jpeg")
, plot = ggplot2::last_plot()
, width = 11
, height = 8.5
, units = "in"
, dpi = "print"
)Adding in bar for remaining area:
constrained_by_scnro_ls_long %>%
dplyr::mutate(
constraint = forcats::fct_recode(constraint,"Available"="Total")
, pct_rdctn_hack = ifelse(constraint=="Available", 1+pct_rdctn, pct_rdctn)
) %>%
dplyr::mutate(
dplyr::across(
c(constraint)
, ~ forcats::fct_rev(.x)
)
) %>%
# plot
ggplot(
mapping = aes(
y = scenario_lab
, x = pct_rdctn_hack
, fill = constraint
, group = constraint
)
) +
geom_col(
width = 0.7, alpha = 0.9
) +
geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
geom_text(
mapping = aes(
label = scales::percent(
ifelse(pct_rdctn<0.145*-1 & constraint!="Available",pct_rdctn*-1,NA)
, accuracy = 1
)
, group = constraint
, fontface = "bold"
)
, position = position_stack(vjust = 0.5)
, color = "black", size = 1.8
) +
geom_text(
mapping = aes(
label = scales::percent(
ifelse(constraint=="Available",pct_rdctn_hack,NA)
, accuracy = 1
)
, group = constraint
, fontface = "bold"
)
, color = "black", size = 2.7
, hjust = -0.1
) +
annotate(
geom = "text"
, x = -0.78
, y = 0
, label = "constrained"
, color = cols_cnstrnt_avail[1]
, fontface = "bold"
, size = 2.5
) +
annotate(
geom = "text"
, x = 0.88
, y = 0
, label = "available"
, color = cols_cnstrnt_avail[2]
, fontface = "bold"
, size = 2.5
) +
scale_fill_manual(values = c(rev(plasma_custom), cols_cnstrnt_avail[2])) +
scale_x_continuous(
limits = c(-1.00,1.02)
, labels = scales::percent_format()
) +
scale_y_discrete(expand = expansion(mult = c(0.8,0.3))) +
facet_wrap(
facets = vars(area_name)
, ncol = 3
, labeller = label_wrap_gen(width = 35, multi_line = TRUE)
, scales = "free_y"
) +
labs(
fill = ""
, x = "Percent Available and Constrained Relative to Total Forest and Shrub Area"
# , x = "Percent of Forest & Shrub Area in Overall Landscape"
, subtitle = ""
, y = ""
) +
theme_light() +
theme(
legend.position = "top"
, legend.direction = "horizontal"
, legend.title = element_text(size=7)
, axis.title.x = element_text(size=10, face = "bold")
, axis.title.y = element_blank()
, axis.text.x = element_blank()
, axis.ticks.x = element_blank()
, panel.grid.major = element_blank()
, panel.grid.minor = element_blank()
, strip.text = element_text(color = "black", face = "bold", size = 9)
) +
guides(fill = guide_legend(nrow = 1))
# export
ggplot2::ggsave(
filename = paste0("../data/plt_cnstrnt_avail.jpeg")
, plot = ggplot2::last_plot()
, width = 8.5
, height = 11
, units = "in"
, dpi = "print"
)constrained_by_scnro_ls_long %>%
dplyr::mutate(
constraint = forcats::fct_recode(constraint,"Available"="Total")
, pct_rdctn_hack = ifelse(constraint=="Available", 1+pct_rdctn, pct_rdctn)
) %>%
dplyr::mutate(
dplyr::across(
c(constraint)
, ~ forcats::fct_rev(.x)
)
) %>%
dplyr::mutate(
dplyr::across(
c(scenario_lab, area_name)
, ~ forcats::fct_rev(.x)
)
) %>%
# plot
ggplot(
mapping = aes(
y = area_name
, x = pct_rdctn_hack
, fill = constraint
, group = constraint
)
) +
geom_col(
width = 0.7, alpha = 0.9
) +
geom_vline(xintercept = 0, color = "black", linetype = "solid", linewidth = 1.1) +
geom_text(
mapping = aes(
label = scales::percent(
ifelse(pct_rdctn<0.145*-1 & constraint!="Available",pct_rdctn*-1,NA)
, accuracy = 1
)
, group = constraint
, fontface = "bold"
)
, position = position_stack(vjust = 0.5)
, color = "black", size = 2
) +
geom_text(
mapping = aes(
label = scales::percent(
ifelse(constraint=="Available",pct_rdctn_hack,NA)
, accuracy = 1
)
, group = constraint
, fontface = "bold"
)
, color = "black", size = 3
, hjust = -0.1
) +
annotate(
geom = "text"
, x = -0.72
, y = (constrained_by_scnro_ls_long %>% dplyr::pull(area_name) %>% unique() %>% length())*1.035
, label = "constrained"
, color = cols_cnstrnt_avail[1]
, fontface = "bold"
, size = 3
) +
annotate(
geom = "text"
, x = 0.82
, y = (constrained_by_scnro_ls_long %>% dplyr::pull(area_name) %>% unique() %>% length())*1.035
, label = "available"
, color = cols_cnstrnt_avail[2]
, fontface = "bold"
, size = 3
) +
scale_fill_manual(values = c(rev(plasma_custom), cols_cnstrnt_avail[2])) +
scale_x_continuous(
limits = c(-1.00,1.02)
, labels = scales::percent_format()
) +
scale_y_discrete(expand = expansion(mult = c(0.02,0.05))) +
facet_grid(cols = vars(scenario_lab)) +
labs(
fill = ""
, x = "Percent Available and Constrained Relative to Total Forest and Shrub Area"
# , x = "Percent of Forest & Shrub Area in Overall Landscape"
, subtitle = ""
, y = ""
) +
theme_light() +
theme(
legend.position = "top"
, legend.direction = "horizontal"
, legend.title = element_text(size=7)
, axis.title.x = element_text(size=10, face = "bold")
, axis.title.y = element_blank()
, axis.text.x = element_blank()
, axis.text.y = element_text(color = "black",size=10, face = "bold")
, axis.ticks.x = element_blank()
, panel.grid.major = element_blank()
, panel.grid.minor = element_blank()
, strip.text = element_text(color = "black", face = "bold", size = 11)
) +
guides(fill = guide_legend(nrow = 1))
4.5 Fireshed-level analysis
Based on model simulations of how much area generally needs to be treated to influence wildfire behavior, we binned the firesheds into three classes of mechanical constraint:
- high (81–100% [i.e., only 0–19% is available for mechanical treatment]): fuels treatment would principally need to rely on fire
- medium (60–80% [i.e., 20–40% is available]): could use a combination of fire and mechanical thinning
- low (60% [i.e., >40% is available]): could effectively influence wildfire behavior with mechanical treatment alone
4.5.1 Distribution of Fireshed PA Constraint
constrained_by_scnro_ls_pa %>%
dplyr::count(area_name,scenario_lab,cnstrnt_class) %>%
dplyr::group_by(area_name,scenario_lab) %>%
dplyr::mutate(
pct = n/sum(n)
, high_pct=max(ifelse(cnstrnt_class=="high constraint",pct,0))
, tot = sum(n)
, scenario_lab = scenario_lab %>% forcats::fct_rev()
) %>%
ggplot() +
geom_col(
mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), fill=cnstrnt_class)
,width = 0.7, alpha=0.8
) +
geom_text(
mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), group=cnstrnt_class
,label = scales::percent(ifelse(pct>=0.12,pct,NA), accuracy = 1)
, fontface = "bold"
)
, position = position_stack(vjust = 0.5)
, color = "black", size = 2.3
) +
# geom_text(
# data = constrained_byftr_huc12_wide %>%
# dplyr::group_by(area_name) %>%
# dplyr::summarise(n=n(), high_n=sum(ifelse(cnstrnt_class=="high constraint",1,0)))
# , mapping = aes(x = -0.05,y=reorder(area_name, -(high_n/n)),label=paste0("n=",scales::comma(n)))
# , size = 3
# , color = "black"
# , hjust = 0.7
# ) +
scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
facet_grid(cols = vars(scenario_lab)) +
scale_x_continuous(labels = scales::percent_format()) +
labs(
fill = "Level of\nConstraint"
, y = ""
, x = "Percent of Fireshed Project Areas"
) +
theme_light() +
theme(
legend.position = "top"
, legend.direction = "horizontal"
, legend.title = element_text(size=7)
, axis.title.x = element_text(size=10, face = "bold")
, axis.title.y = element_blank()
, axis.text.x = element_blank()
, axis.text.y = element_text(color = "black",size=10, face = "bold")
, axis.ticks.x = element_blank()
, strip.text = element_text(color = "black", face = "bold", size = 11)
) +
guides(
fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
)
4.5.2 Distribution of Fireshed Area
# plot
ggplot(data = constrained_by_scnro_ls_pa %>% dplyr::filter(scenario_id==1)
, mapping = aes(x = feature_area_ha, y = reorder(area_name,desc(area_name) ))
) +
geom_vline(
xintercept =
constrained_by_scnro_ls_pa %>%
dplyr::filter(scenario_id==1) %>%
dplyr::pull(feature_area_ha) %>%
median()
, linetype="dashed", color="gray66"
) +
geom_violin(aes(fill = state)) +
geom_boxplot(width = 0.15) +
# geom_point(size = 0.7, color = cols_obj_r[2], alpha = 0.2) +
geom_text(
data = constrained_by_scnro_ls_pa %>%
dplyr::filter(scenario_id==1) %>%
dplyr::group_by(area_name) %>%
dplyr::summarise(n=n(), max_area=max(feature_area_ha))
, mapping = aes(
x = (
constrained_by_scnro_ls_pa %>%
dplyr::filter(scenario_id==1) %>%
dplyr::pull(feature_area_ha) %>%
min()
)*.9
,y=reorder(area_name,desc(area_name) )
,label=paste0("n=",scales::comma(n))
)
, size = 3
, color = "black"
, hjust = 0.7
) +
scale_fill_viridis_d(option = "viridis", alpha = 0.8) +
scale_x_continuous(labels = scales::comma, breaks = scales::extended_breaks(n=7)) +
labs(
fill = ""
, y = ""
, x = "Fireshed Area (ha)"
) +
theme_light() +
theme(
legend.position = "none" # c(0.9, 0.9)
, axis.title = element_text(size=9)
, axis.text.x = element_text(size=7)
)
4.5.3 Map of Firesheds by Level of Constraint
# polish for mapping
map_firesheds <- fireshed_proj_area %>%
dplyr::select(pa_id, geometry) %>%
dplyr::mutate(pa_id=as.character(pa_id)) %>%
dplyr::inner_join(
constrained_by_scnro_ls_pa
, by = dplyr::join_by("pa_id")
, multiple = "all"
) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("pct_")
, ~ scales::percent(as.numeric(.x), accuracy = 0.1)
)
, dplyr::across(
tidyselect::ends_with("_ha")
, ~ scales::comma(as.numeric(.x), accuracy = 1)
)
, scenario_lab = scenario_lab %>% forcats::fct_rev()
) %>%
dplyr::mutate(
Priority.Area.Name = area_name
, Level.of.Constraint = cnstrnt_class
, Pct.Treatable.Class = rmn_cnstrnt_class
, Fireshed.ProjArea.ID = pa_id
, Fireshed.ProjArea.Area.Hectares = feature_area_ha
, Fireshed.WCS.Status = ifelse(is.na(fireshed_crisis_strategy),"Not High Risk",fireshed_crisis_strategy)
, ForestShrub.Area.Hectares = covertype_area_ha
, Pct.Rdctn.Protected = pct_rdctn1_protected
, Pct.Rdctn.Slope = pct_rdctn2_slope
, Pct.Rdctn.Roads = pct_rdctn3_roads
, Pct.Rdctn.Riparian = pct_rdctn4_riparian
, Pct.Rdctn.Administrative = pct_rdctn5_administrative
, Treatable.Area.Hectares = rmn5_administrative_area_ha
, Pct.Treatable.Area = pct_rmn5_administrative
)
# filter for mapview
mapview_fs_temp <- map_firesheds %>%
dplyr::filter(
scenario_id == 2
)
# ## filter for high priority only
# hi_map_firesheds <- map_firesheds %>%
# dplyr::filter(!is.na(fireshed_crisis_strategy))
# filter for mapview
mapview_fs_sc1_temp <- map_firesheds %>%
dplyr::filter(
scenario_id == 1
)
mapview_fs_sc3_temp <- map_firesheds %>%
dplyr::filter(
scenario_id == 3
)
mapview_fs_sc3_hi_temp <- map_firesheds %>%
dplyr::filter(
scenario_id == 3
& fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
)# mapview
mapview::mapviewOptions(homebutton = FALSE, basemaps = c("OpenStreetMap","Esri.WorldImagery"))
mapview::mapview(wf_landscapes
, color = "black"
, lwd = 1
, alpha.regions = 0
, label = FALSE
, legend = FALSE
, popup = FALSE
, layer.name = "WCS Landscapes"
, hide = TRUE
) +
mapview::mapview(
mapview_fs_temp
, zcol = "cnstrnt_class"
, col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
, alpha.regions = 0.6
, lwd = 0.2
, label = FALSE
, legend = FALSE
, layer.name = "Scenario 2 Constraints"
, hide = FALSE
, popup = leafpop::popupTable(
mapview_fs_temp
, zcol = c(
"Priority.Area.Name"
, "Level.of.Constraint"
, "Pct.Treatable.Class"
, "Fireshed.ProjArea.ID"
, "Fireshed.ProjArea.Area.Hectares"
, "Fireshed.WCS.Status"
, "ForestShrub.Area.Hectares"
, "Pct.Rdctn.Protected"
, "Pct.Rdctn.Slope"
, "Pct.Rdctn.Roads"
, "Pct.Rdctn.Riparian"
, "Pct.Rdctn.Administrative"
, "Treatable.Area.Hectares"
, "Pct.Treatable.Area"
)
, row.numbers = FALSE
, feature.id = FALSE
)
) +
mapview::mapview(
mapview_fs_sc1_temp
, zcol = "cnstrnt_class"
, col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
, alpha.regions = 0.6
, lwd = 0.2
, label = FALSE
, legend = FALSE
, layer.name = "Scenario 1 Constraints"
, hide = TRUE
, popup = leafpop::popupTable(
mapview_fs_sc1_temp
, zcol = c(
"Priority.Area.Name"
, "Level.of.Constraint"
, "Pct.Treatable.Class"
, "Fireshed.ProjArea.ID"
, "Fireshed.ProjArea.Area.Hectares"
, "Fireshed.WCS.Status"
, "ForestShrub.Area.Hectares"
, "Pct.Rdctn.Protected"
, "Pct.Rdctn.Slope"
, "Pct.Rdctn.Roads"
, "Pct.Rdctn.Riparian"
, "Pct.Rdctn.Administrative"
, "Treatable.Area.Hectares"
, "Pct.Treatable.Area"
)
, row.numbers = FALSE
, feature.id = FALSE
)
) +
mapview::mapview(
mapview_fs_sc3_temp
, zcol = "cnstrnt_class"
, col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
, alpha.regions = 0.6
, lwd = 0.2
, label = FALSE
, legend = FALSE
, layer.name = "Scenario 3 Constraints"
, hide = TRUE
, popup = leafpop::popupTable(
mapview_fs_sc3_temp
, zcol = c(
"Priority.Area.Name"
, "Level.of.Constraint"
, "Pct.Treatable.Class"
, "Fireshed.ProjArea.ID"
, "Fireshed.ProjArea.Area.Hectares"
, "Fireshed.WCS.Status"
, "ForestShrub.Area.Hectares"
, "Pct.Rdctn.Protected"
, "Pct.Rdctn.Slope"
, "Pct.Rdctn.Roads"
, "Pct.Rdctn.Riparian"
, "Pct.Rdctn.Administrative"
, "Treatable.Area.Hectares"
, "Pct.Treatable.Area"
)
, row.numbers = FALSE
, feature.id = FALSE
)
) +
mapview::mapview(
mapview_fs_sc3_hi_temp
, zcol = "cnstrnt_class"
, col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
, alpha.regions = 0.6
, lwd = 0.2
, label = FALSE
, legend = FALSE
, layer.name = "Scenario 3 Constraints (High Only)"
, hide = TRUE
, popup = leafpop::popupTable(
mapview_fs_sc3_hi_temp
, zcol = c(
"Priority.Area.Name"
, "Level.of.Constraint"
, "Pct.Treatable.Class"
, "Fireshed.ProjArea.ID"
, "Fireshed.ProjArea.Area.Hectares"
, "Fireshed.WCS.Status"
, "ForestShrub.Area.Hectares"
, "Pct.Rdctn.Protected"
, "Pct.Rdctn.Slope"
, "Pct.Rdctn.Roads"
, "Pct.Rdctn.Riparian"
, "Pct.Rdctn.Administrative"
, "Treatable.Area.Hectares"
, "Pct.Treatable.Area"
)
, row.numbers = FALSE
, feature.id = FALSE
)
)